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I. 



INTRODUCTION AND SUMMARY 



A. INTRODUCTION 

The failure rate function, or hazard function ( hazard 
for short) may be described as the conditional probability 
of an equipment's failing at operating age t, having sur- 
vived to that age. The reliabilities of a variety of elec- 
tronic and mechanical items are conveniently and naturally 
described in terms of the appropriate hazard fxinction, and 
so is the longevity of human beings. The term force of 
mortality replaces hazard in the latter context. 

This paper is devoted to a study of several simple 
analytical representations for hazard functions. These 
representations are in turn based upon representations of 
random variables having certain required properties, in 
terms of others having familiar distributions — in particular 
the exponential. Similar ideas are due to Tukey [1] , and 
recently have been examined by Parzen [2] . The hazard 
representations proposed are quite expeditiously used in 
simulation studies, e.g. of system reliability or avail- 
ability in terms of component lifetimes. They may also be 
used in data analysis studies, in order to parsimoniously 
describe data sets in terms of perturbations of convenient 
and familiar standard distributions. Their use in data 
analysis and simulation is also described in Gaver, Laven- 
berg, and Price [3], and in Gaver and Chu [4]. 
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B. SYSTEM FAILURE PATTERNS 



It is plausible to think that the time series of fail- 
ures in a system may involve these stages. 

Early failures . There may be a relatively large number 
of failures soon after a system is introduced because of 
design defects, production errors, or errors stemming from 
maintenance personnel inexperience. This situation is 
characterized by a hazard function that is initially large, 
but that decreases with time. "Infant mortality" is in 
evidence . 

Random Fai lures . Following the early failure period 
there may be a period during which failures occur at an 
essentially constant rate for a rather prolonged time. 
During this period the hazard function is nearly constant, 
so the times between failures are close to being exponen- 
tially distributed. The effect of age or wearout is not 
yet apparent. 

Wearout Failures . Eventually following the period 
during which a constant hazard is evident there is likely 
to be a period of ever-increasing failure rate caused by 
wearout of system components. 

A graphical representation of a hazard function that 
exhibits the behavior described is given below. Note that 
it has the legendary "bath tub" shape. 
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Fig. 1. Bath Tub Shape Curve 



Some comments on the above follow: 

The term "failure" may refer to an event that is 
analogous to human death, after which the entire system 
is replaced. On the other hand repair or component replace- 
ment may occur after failure: the system is only repaired, 

not entirely replaced. In the fomer case, a hazard func- 
tion of the kind depicted in Figure 1 applies to each system 
event ("death"); when the system is installed (or is born), 
that hazard operates starting from scratch at t = 0 until 
system failure (hviman death, for instance) , after which a 
similar hazard goes into effect, starting once again from 
zero. In the latter case, in which repair of a component 
occurs, a hazard function like that of Figure 1 applies at 
t = 0, but after the first event ("failure") at t^ a re- 
pair action is accomplished. The same hazard operates for 
t >_ t^ until the next event at > tj^, and so on. Inter- 

mediate situations may be envisioned, in which after event 
n at t^ the hazard governing system failure n+1 starts 

at t - T,0<T <t. 

n n n n 

Although there is reason to assume that hazards some- 
what like that of Figure 1 occur in general for systems , 
the possibility exists that the system hazard is "bumpy" 
because wearout failures of components or subsystems may 
well occur at intermediate times. 

If the theory is applied to systems with little or no 
wearout propensity, as should be the case when dealing with 
computer software modules, then the hazard function may 
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well exhibit the initial falloff of Figure 1 but not the 
rise at later times. In fact, a constant decline as bugs 
are fo\md and removed could be (optimistically) anticipated 
for software. The right-hand side of the bath tub vanishes, 
and the picture is that of a ski slope. 
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II. ANALYTICAL HAZARD REPRESENTATIONS 



A. MODELS FOR THE HAZARD FUNCTION 

In this section mathematical models are presented for 
the failure rate or hazard function. Recall that the hazard 
may be defined as follows. 

Definition . Suppose that the time to failure, X, is a 
random variable with distribution function F(x), where 
F(0) = 0; the latter possesses the density function f(x), 
f(x) = dF/dx, such that for any positive x. 



X 

F(x) = / f(y) dy . (2.1) 

0 



Then the hazard function, or failure rate at age 
given by 



h (x) 



f (X) 

1 - f(x) 



X, is 



( 2 . 2 ) 



The interpretation of h(x)dx is that it is the con- 
ditional probability of failure in the interval (x, x + dx) , 
given that there has been no failure up to age x. 

Express the hazard as 



h(x) 
h (x) dx 



dF/dx 
1 - F(x) 

dF 

1 - F(x) 



- d{ log [1 - F (x) ] } , 



it then follows after integration that 
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(2.3) 



X 

F(x) = 1 - exp(- / h(y)dy] . 

a 

Thus if the hazard is specified, so is the distribution 
function, and conversely. 

Note that if 



h(x) = X > 0 , 


0 £ X 




F(x) = 1 - e , 


0 < X , 


(2.4) 



so a constant hazard fxanction implies the exponential dis- 
tribution of the random variable X, and conversely. 

Obviously a constant hazard representation does not 
describe the bath tub hazard shape of Figure 1, nor does it 
represent a situation in which hazards decline, possibly 
because design defects or "bugs” are occasionally removed. 
Here are two hazard representations likely to be useful for 
such purposes. 

1 . A Bath Tub Model 

Define the random variable Z in terms of X, X 
being exponentially distributed with mean X as follows: 

Z = G(X) = XL(X)R(X) 

or (2.5) 

= X(J>(X) , (j)(X) = L(X)R(X) 

where 

a) L(x) is concave in x, L(0) < 1, L(«) = 1, 

b) R(x) is convex in x, R(0) = 1, R(0) > R(<») . 
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Then the hazard of Z may be made to exhibit a bath t;ab 
shape, as in Figure 1, by proper choice of the functions L 
and R. 



Example . Suppose 



L(x) = 



ax 

1 + ax 



a > 0 , 0 < X 



( 2 . 6 ) 



R(x) = 



1 + 3x 



3 > 0 



Clearly, 



z = xL(x)R(x) = X 



ax 



ax 



X 



1 + ax l + 3x 1 + ax l+3x 



is a monotonically increasing function of x. Furthermore, 
choose a large (e.g. a = 10) and 3 small (e.g. 3 = 10 ^) . 
Then it is intuitively clear that (i) small x- values 
transform into even smaller z-values, e.g. x = 1 corre- 
sponds to z = 0.91 and x = 2 corresponds to z = 1.90, 
but (ii) this effect dwindles as x increases, so x = 10 
corresponds to z = 9.8 and x=50to z=47.5 and the 
z-values closely resemble the x's percentage-wise, but 
(iii) as X increases still further the z's do not follow 
suit: X = 10^ corresponds to z = 500. This suggests 

that if X is a value assumed by X, that Z shares the 
properties of X in mid-range, i.e. for intermediate 
x-values, but differs from X by having a disproportionate 
probability of assuming small values (near zero) , or large 
values (near, but less than, 1/3) • Thus the hazard of 
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Z will appear to be a "bath tubbed" version of X, 
particularly if X is exponential. 

We focus attention on the representation (2.6) in what 
follows, mainly for analytical and computational convenience. 
Of course there are many other possibilities, such as 

L(x) = 1 - e"“^ 

-Bx 

R(x) = e ; 

these latter may be adjusted to provide sharper-edged tubs 
than can (2.6), but iteration of (2.6) may be induced to 
accomplish the same purpose. 

2 . A Decreasing Failure Rate Model 

Define the random variable W in terms of X, X 
again being exponential with parameter 1 

W = XT(X) (2.8) 

where T(x) is an increasing function of x, L(0) = 1. 

Then the hazard may be made to exhibit a decreasing 
behavior. 



Example . Suppose 



T(x) = 1 + cx , 



c > 0 , 0 < X . 



(2.9) 



Then 



2 = x(l + cx) 



( 2 . 10 ) 



is monotonic, and small x-values lead to comparable z-values 
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(especially when c is small) , but larger x-values are 
"amplified” by 1 + cx to yield increasingly large z-values. 

Attention will be focused upon (2.9), although other 
possibilities exist that accomplish the same purpose, namely 
that of lengthening the right tail of the distribution of 
X (simulating outliers, for instance) while leaving the 
body of the distribution virtually unchanged. 

B. MATHEMATICAL PROPERTIES OF THE "BATH TUB" HAZARD MODEL 

Various analytical properties of the previously described 
models will now be recorded. These provide useful insights 
into the behavior of the random variables Z and the under- 
lying (generating) variables X. 

1. Monotonicity; Quantiles 

It is convenient to focus on monotonic increasing 
transformations, i.e. if 



z = G(z)=Z(|)(z) (2.10) 

then in order that the above function be monotonically in- 
creasing, dz/dx > 0. Observe that logarithmic differen- 
tiation of (2.10) provides 



dz 

z 



dx (p‘ (x) 

X (t> (x) 



dx 



( 2 . 11 ) 



and thus dz/dx >0 if and only if 



1 . 4>' (X) 
X 4> (x) 



> 0 



( 2 . 12 ) 
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Alternatively, the condition is, in terms of L(x) and 
R(x) , 



1 , L' (X) , R' (X) 

X L (x) R(x) 



> 0 



(2.13) 



It is easily seen that the important example (2.6), 



<P (x) = 



ax 



1 + ax 1 + 6x ' 



yields a monotonic relationship between z and x. The 
fact that this transfomiation can be easily and explicitly 
inverted (solved for x in terms of z) will be exploited 
subsequently. 

Of course if z(x) is monotonically increasing then so 
is x(z) , the inverse function. The events (Z £ z) and 
(X £ x(z)) are equivalent, and so 

P{Z £ z} = P{X £ x(z)} , (2.14) 

from which it follows that if x = x(p) is the p*100% 

P 

quantile of X, i.e. 

p{X <_ x(p) } = p , (2.15) 

then 

P{Z <_ z (p) } = P{Z <_ z(x(p))} = p (2.16) 

and so z(p), the p*100% quantile of Z is simply obtained 
from 

z(p) = x(p) (|)(x(p)) = x(p) L(x(p)) R(x(p)) (2.17) 

In other words we very easily translate from (points on) 
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the inverse distribution of X to the inverse distribution 



of Z . Explicit representation of the distribution of Z 
is however, not often easily possible. 



2 . Hazard and Density Function Relationships 

In order to investigate the relationship between 
the hazards of Z and X, begin by writing 

x(p) 

p = F„(x(p)) = 1 - exp[- / h^(u)du] (2.18) 

X Q X 



or 

x(p) 

/ h (u)du = - Jln(l-p) 

0 ^ 



Now differentiate with respect to p to find 



SgE) 



or 



= difpr • = 



(2.19) 



( 2 . 20 ) 



here h and f are the hazard and density functions of 

X X 

the r.v. X. The relationship (2.20) holds for any distri- 
bution, of course. 

Differentiation of (2.5) reveals the connection between 

h and h . From (2.11) 
z X 



dz (p) 
dp 



= z (p) 



1 , 0' (x(p) ) 

x(p) <p (x(p) ) 



dx(p) 

dp 



( 2 . 21 ) 



[<j)(x(p)) + x(p) <{)'(x(p))] 
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n 



From (4.11), applied now to the z-hazard, there results 



h (zTpT) ° hTxtp)) l**^*?)’ ♦'(X(P))1 , (2.22) 

so 

h (z(p)) = h (x(p)) XT — 7 — TT — ; 7 —r — m — ? — rr- (2.23) 

z ^ X ^ (}) (x (p) ) + x(p) (x(p) ) ' 

Multiplication of both sides by 1-p then shows, in view 
of (4.11), that the density functions are similarly related 

f2<2(p)) = fj;(x(p)) 4 ,(x'(p)) + x(pt (x(p))- • 



Example . 

X is exponential (A ) . Then 



^^(^(p)) (}) (x(p) ) + x(p) (})' (x(p) ) 



(2.24) 



Now use the specific 0(x) of (2.6) 



.L / \ _ cix 

■ 1 + ax ‘ 1 + 6x ' 



or, in terms of logarithms. 



in 4>(x) = in ax - in(l + ax) - in(l + 3x) , 



so 



<})' (x) 
4> (x) 



_ 1 _ g _ 1 - ggx 25) 

X 1 + ax 1 + 6x x(l + ax) (1 + 3x) ' 



and 



(f) (X) + X(j)' (x) = (p (x) 



2 + (g + 3 ) X 
(1 + ax) (1 + 3x) 



(2.26) 



finally 
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(2.27) 



h (2(p)) 



(j)(l + gx(p))^ (1 + 3x(pn^ 
ax(p) [2 + (a + B) x(p) ] 



Although this expression is not quite explicit, qualitative 
properties of h^ can be deduced from it. 

If p ->■ 0 , x(p) = " iind-p) + 0, and hence 



or 



h (z(p)) ~ -s — i- 7 — p, 
z' ' 2a x(p) ' 



lim x(p) h^(z(p) ) = j- 

p 0 



Since for p ^ 0, 



z (p) ~ ax (p) 



and hence 



x(p) ~ [z(p)/a] 



1/2 



there results 



h (z(p)) ~ ^ 1 / 2 ' ~ 

^ 2[z(p) /^ 



or 



lim /z (p) h (z(p)) = 

p ^ 0 2/a' 



(2.28) 

(2.29) 



(2. 30) 



(2.31) 



This shows that h (z) f « as z 0 , creating the left- 

z 

hand end of the bath tub of Figure 1. 

If p 1, x(p) f «, and z(p) f 1/8 so 



h (z (p) ) 

2 



~ A 



a8^ x^ (p) 

Xa + 8) 



or 
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(2.32) 



lim 

p oo 




For p -»■ 1 




(2.33) 



so 




(2.34) 



and thus 



h (z (p) ) ~ X 




(1 - 32 )^ 



1 



(2.35) 



Once again it appears that the hazard rises rapidly, this 



The bath tub effect is presiomably achieved by choosing a 
large and 3 small. Let a 0 and 3 0 independently 

in (2.36); it is clear that the limiting value of the 
hazard is X . This indicates that the hazard is (approxi- 
mately) X for middling values of z. 

3 . An Explicit Formula for a Hazard 

The expression (2.6) leads to the relationship 



time as x(p) t and z(p) t 3 the other end of the 



bath tub is thus fashioned. 

If p = 1 - e then x(p) = X Then 



\{z{l - e~h) 



= [1 + g/X]^ [1 + 3/X]^ (2.36) 

a/X [2 + (a + 3)/X] 



" (1 + ax(pJT (1 + 3x(pyy 




(2.37) 
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and the latter may be explicitly inverted by solving a 
quadratic equation. The result is 



x(p) = (o^'*'S)z(p) + v^(g+S) 



2a(l - 



z (p) + 4z (p) 
Bz (p) ) 



.“.(1- Bz (£)).., (2.38) 



Now a direct differentiation of this expression and invo- 
cation of (2.20) produces the expression 



h^(z) 



jQTrr-BzT 



6{(a+6)z+ '/(a+6)^z^ + 4az(l-6z) , , , ^\ 

+ (a + 6) 



{ (g-B) ^z + 2a}y(g+g) ^z^ + 4az (l-gz) 
(a + B) z + 4az (1 - 6z) 



(2.39) 



This form, while explicit, provides no particularly useful 

insights; the bath tub end shapes already noted in (2.31) 

and (2.35) can be deduced directly from (2.39). 

Some graphical plots of h are presented below. They 

z 

illustrate the behavior of the present hazard representation 
in a more understandable fashion than does the formula 
itself. 

4 . An Explicit Formula for the Failure Time Distribution 
Because x and z are monotonically related 
through (2.37) we have 
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p{z < z} 



( 



+ y / - . « ' 2 2 



_ r. ) V / _ (ci+6) z + V (a+6) z + 4az (1-Bz) 

_ P |X < X 2a U"- 6z) 



( 



= 1 “ exp \ -X 



V ' ' 2 _2 



(a+6)z + V (a+B) z + 4az(l-Bz) 
2a(l - Bz) 



(2.40) 



Again the explicit formula seems unproductive of insights. 



C. MATHEMATICAL PROPERTIES OF THE DECREASING FAILURE 
RATE MODEL 

1 . The Hazard Behavior 

The expression (2.23) can be applied to deduce the 
hazard function of the representation (2.8), advertised to 
produce a decreasing failure rate. There we specified 



({> (x) = T (x) = 1 + cx , 



(2.41) 



and thus, from (2.2 3) 



hw(w(p)) (1 + cx (p) ) + x(p) c 1 + 2cx(p) 



(2.42) 



Qualitative properties follow easily. 

If p 0 , x(p) 4 0, and 

h (w(p) ) ~ X 
w 

Thus the hazard is approximately X for small z. 
If p 1, x(p)' f and 



(2.43) 



w(p) ~ c [x(p) ] 



so 
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h ,(w(p) ) 



w 



(2.45) 



2 /cw ( p ) 

which clearly decreases, as claimed. It may be inferred 
that the distribution of W appears nearly exponential , 
but possesses an extraordinarily long right tail — possibly 
the result of outliers. 



2 . Explicit Formulas for the Hazard and the 
Distribution Function 

Direct solution of the quadratic equation 



w = x(l + cx) = X + cx‘‘ 



presents 



, . /I + 4cw(p) - 1 
X(p) = 



which, when differentiated, leads to 



(2.46) 



h (w) = 



w 



/I + 4cw 



h. (x) 



X 



A 



/I + 4cw 



(2.47) 



The distribution function is 

p{w 



X 1 « ) X /I + 4cw - 1 ( 

<w)=pjx<x = ^ 



= 1 - exp 






/I + 4cw - 1 1 



2c 



(2.48) 



This distribution bears a close family resemblance to tlie 
Weibull distribution 1 - F(w) = exp{-kv^}, especially for 
large (right tail) values of w. 



D. AN ALTERNATIVE "BATHTUB" HAZARD REPRESENTATION 

The simple parametric model (2.6) leading to a bathtub- 
shaped hazard is by no means the only possibility. We next 
describe another simple approach. It is that of defining 
a hazard function having an appropriate shape , and then 
deducing the corresponding distribution function, and a 
procedure for sampling from it, rather than proceeding in 
reverse order, as before. 

Let the hazard be of the form 



h(z) = g(z) +•♦•+ k(z) , (2.49) 

where g(z) >0 is a decreasing function of z such that 
lim^ ^ g(z) =0, and k(z) is an increasing function of z, 
such that (preferably) k(0) =0 and k(°o) = ". Such a 
function can yield a bathtubbed hazard. 



Example . 

+ BZ + X 



(2.50) 



A, B, a, A all positive. 

Clearly, (2.50) has a generally "bath tub-like" appear- 
ance, since 



if 



while for 



h(z) >0. 



h' (z) — ^ + B (2.51) 

(z + a) ^ 

- - + B < 0, then h' (0) < 0 (2.52) 

a 

z > Zq = g - a (2.53) 
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Detailed behavior is adjustable by choice of the 
parameters. Now the distribution function of time to 
failure, Z, is obtained from (2.50) : 



P{Z > z} = 



( " 

exp<- / h(x)dx 
( 0 






exp } - I A in (1 + ^ + Xz 



I 



[ 



2v . B _2 
a' 2 



] 



A 

/ ot X / B 2^ -A.Z 



= F^(z) F 2 (z) F^Cz) , 



(2.54) 



where 



■ ( rfr)* 

F 2 (z) = exp f 
F3(z) = e~^^. 



(2.55) 



All of the above are recognized as being the complements of 
distribution functions. In effect, the distribution of Z 
is that of the minimum of three independent random variables 

P{Z > z} = P{X^ > z}*P{X 2 > z}*P{X 3 > z}, (2.56) 
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having the distribution = 1 - (i = 1,2,3). 

This fact leads directly to an easy procedure for simulation 
of Z by simply obtaining the smallest from among the 
realization of X^, and X^. The advantage of the pre- 

vious method, based on (2.6) for instance, is that only one 
realization — that of an exponential in that specific example 
— leads to the realization of Z. This is not only computa- 
tionally attractive, but seems to facilitate the application 
of such Monte Carlo variance reduction techniques as control 
and antithetic variables, cf., Hammersley and Handscomb [5]. 
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III. OBTAINING SPECIFIED HAZARD BEHAVIOR 



BY SIMPLE SAMPLING 



The development of the last section illustrates one 
manner in which hazard behavior may be conveniently repre- 
sented and simulated. We now show how such behavior may 
alternatively be obtained by simple simulation, i.e. from 
one realization of a basic (possibly exponential) random 
variable . 

Refer to (2.5), in which 



Z = G(X) (3.1) 

and, if G(*) is monotonically increasing, 

z(p) = G(x(p)) , ^ (3.2) 

z(p) and x(p) being the p*100% percentiles of Z and X, 
respectively. Then the counterpart to (2.23) that 'results 
from differentiation of (3.2) is the expression 

h^(2(p)) = h^(x(p)) = h^(x(p)) 

Consequently, if one specifies h (z) as a suitable func- 
tion of the "time" z, and specifies the distribution of 
the stochastic variable X — and hence its hazard, h^ — there 
results a differential equation for z(x) = G(x) ; 



h^(2) 



dz 

dx 



= h^(x) 



(3.4) 



integration then provides the desired transformation, G. In 
other words, we seek z(x) satisfying 



28 



r 



(3.5) 



2 X 

/ h (u)du = / h (v)dv 
0 ^ 0 ^ 



which can sometimes be carried out in a useful closed form. 

Example 3.1 . Refer to the example of Section II, wherein 
h^ is given by the expression (2.50) and we assume that X 
is exponential, so h^ is constant. Then in order to deter- 
mine G(x) = z(x), solve the equation 



or 



/ 

0 



u + a 



+ Bu + A. 



du = 



X 

/ dv 
0 



A Jln(l + -) + I 
a 2 



+ Xz = X 



(3.6) 



Closed-form solution of this expression for z in terms of 
X is of course impossible. One possible approach is purely 
numerical; find an approximate solution, Zq(x), e.g. the 
appropriate solution of the quadratic 

I z^ + Xz = X (3.7) 

and then correct the result by a few Newton- Raphson itera- 
tions. In other words, put 



2 ^ (x) 



- X + \/x^ + 2Bx 
B 



(3.8) 



now apply Newton to obtain an improved solution 

A in (1 + 2^/a) 

22 (x) = - A/( 2 j_ + a) + B 2 ^ + X 



(3.9) 
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which will b 0 feasible if 0 < Z 2 « The process can be 
iterated (the numerator will change after the first iter- 
ation) . If one wishes to use this model it may actually be 
desirable to start by solving 

I + (A + §)z - X = 0 (3.10) 

^ ot 

for z^, in which case the numerator will not be as shown in 
(3.9); convergence may be more rapid. 



Example 3.2 . Change the hazard representation of the previous 
example as follows; let 



h_(u) = 



(v + a) 



+ Bu + A; (A, a, B, A > 0) (3.11) 



then 



/ h (u)du = ^ ^ + !• z^ + Az (3.12) 

n z a(z + a) 2 



Now it is necessary to solve 



Az 



a ( z + a) 



, B 2 ^ , 

+ j z + Az = 



(3.13) 



i . e . , the cubic 



^ z^ + [a y + A]z^ + [— + aA - x] z - ax = 0 , (3.14) 

which can be carried out, at least formally, in closed form. 
Once again an iterative solution that begins by dropping the 
cubic term, solving the resulting quadratic for z^^(x) , and 
then continuing along the Newton-Raphson road may be^ 
successful. Further investigations of these ideas should 
be conducted. 
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IV. COMPUTER SIMULATION AND ESTIMATION PROCEDURE 



A. SIMULATION AND NUMERICAL RESULTS 

In previous chapters an analytical model was described 
for the failure rate function; useful formulas were also 
derived from the model (2.39) and (2.40). Before the model 
is used in realistic situations, it will be convenient to 
build a computer simulation model for model validation. 

1 . Algorithm 

First a very basic simulation model was built for 
determining the general shape of the failure rate function 
associated with parameters a, 3 and A. In the simulation 
model, a was selected to be 1.0 and 2.0, 3 was selected to 
be 0.05, 0.01, 0.005, 0.001 and X was selected to be 0.1. 
These values were picked arbitrarily, it is only stipulated 
that a is always greater than 3- The general algorithm 
of the simulation model is shown in Fig. 2. 

The hazard function is calculated according to the 
model (2.37) and the system logic function (2.40). The 
results were shown in Fig. 3 and Fig. 4. 

2 . Some Comments 

These simulation results show that : 
a. Parameter a is effective when z values relatively 
have small values. That is, it influences the early 
failure period. 
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FUNCTION 



PLOTTING Z vs FAILURE 
RATE FUNCTION 




Fig. 2. Basic Algorithm. 
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Fig. 3. The Shape of IJazard Function with Different Values 
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Fig. 4. The Hazard Function Over Its Range 



b . 



c . 



d. 



Parameter 3 is effective on the relatively bigger z 
values. That is, it describes wearout failures. 
Parameter X has little effect on the shape of the 
curve; it is a scale factor. 

The last important observation is that the z values 
are limited by the parameter 3 such that; 



" i F 



When z equals 1/3, h (z) goes to infinity, 



B. EMPIRICAL DENSITY FUNCTION AND ESTIMATION PROCEDURE 

In this section, an estimation procedure for parameters 
a, 3, X is defined and the procedure of finding the 
empirical density function is described. 



1. Empirical Density Function 



Suppose that N^, N, A and Nq are defined as 



follows : 



N^ = the frequency of data points for each time interval 
between a^ and b^ for i equals 1,2,3, ...,k. 

N = the total number of failures where. 



N = y N. 

iil ^ 

A = the length of each interval 



A = b . - a . 
1 1 
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Nq = the total number of survivors by the beginning of 
each interval which can be shown as follows: 

N if i = 1 

i-1 

N - I N. if i = 2,3, . . . ,k 

j=l => 

ri=l 

where 2,j = ]_ is the number of failures before 

interval i, or N - V. ^ N. is the number of 

^3=1 3 

failures after interval i. 

If h^(z) is defined as the density function of 
hazard, then the relationship between the frequences of the 
failure data and density function of hazard can be given 
as follows ; 

b. 

1 

~ Nq • / h^(z)dz (4.1) 

a . 

1 

where i = l,2,3,...,k and a^ = (i-1) A, b^ = iA. At this 
point some approximation can be made in expression (4.1), 

such that 

b. 

1 

/ h (z)dz ~ (b. - a.)*h (i(b. - a.)) (4.2) 

^i 

for small values of A = b-a. Then this approximation (4.2) 
is substituted in the expression (4.1), it turns out as 
follows : 

N. = A*N,. h (iA) (4.3) 

1 0 z 

or 
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/ 



1 ^ f f * • n f 



(4.4) 



h^(iA) 



N. 



1 




k 



for the general case, it will be: 



h^(iA) 



N. 

1 

N 



N. 

1 



i-1 

V (N- I N.) 
j = l ^ 



for 



•A 

for 



i 



i 



1 



2,3 



r • • • 



(4.5) 



The expression (4.5) will be used in the calculation of 
an empirical hazard function and also used in the proposed 
estimation procediare. 



2 . Estimation Procedure for Model Parameters 

Two approaches can be used for this problem. The 
first idea is to use the relationship between a sample pth 
percentile and the related probability of the pth percentile. 
A subsequent idea is to approach the problem as a nonlinear 
least square estimation problem for a, 3, X: pick a, 3 
X so that the values obtained minimize the sum of squared 
errors in an objective function. 



a. 3-Percentile Approach 

Suppose F„(x) is the ciomulative probability 
function of the exponential distribution. The pth percentile 
x(p) is equal to the value such that: 

p = Fj^(x(p)) = 1 - (4.6) 
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Then 



Jln(l-p) 



(4.7) 



x(p) = f"^(p) = - j 
where A > 0 . 

If z (p) is defined to be the pth percentile 
in the bathtub model (2.37) then z (p) and x(p) have a 
definite relationship with expression such that; 

= ax (p) 

^ (1 + ax(p)) (1 + 3x(p)) 

by substituting (4.7), 



1 

A 



a{- Zn(l-p) }‘ 



z (p) = 



|l + a [- j iln(l-p) J j j 1 + a [" - j iln(l-p)]j 



(4.8) 



or 



z (p) 



(g/A^)£^(p) 

[l |e(p)}{ 1 + |.(P)| 



(4.9) 



where £ (p) = - in(l-p). 

In the last equation (4.9) , there are three un- 
known variables which are a, 3, A. If three independent 
equations associated with expression (4.9) are defined, 
clearly, they can be solved for the unknown parameters. 
Actually, this can be easily done, for using different values 
of percentiles such that: 
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(a/X^) £^(P^) 

(a/x2)£^(P2) 

{l MP2>} {l + f MPj)} 

‘''3' ■ jlT’iMPj)} {1 + I C(P3)} 

where z(p^) and e(Pj^) are the known values associated 
with P^. Actually, to get the percentile values, two kinds 
of approaches can be made. First they can be computed 
directly from data using the simple statistics method, such 
as 2(j_) (approximately) z(i/(n+l)). Second they can 

be computed from the empirical density function. Generally 
the choice depends on the form of the data that are avail- 
able . 

To decide for the effectiveness of this type of 
estimation, another basic computer simulation is made by 
modifying the first computer simulation algorithm. This 
algorithm is shown in Fig. 5. 

The simulation was rvin four hundred times and 
in each rxm the sample size was assumed to be n = 50. 
Initially the parameters a, B and X were taken to be, 
respectively, 1.0, 0.05 and 0.1. The percentiles P^/ 
and P^ were used in each run such that 0.1, 0.5 and 0.9 
respectively. The estimation results are shown in Table I. 



z(Pj^) = 



2(P2> = 
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START 




Fig . 5 . Simulation algoritm for estimation procedure 
and emprical density function. 
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As can be seen, the estimated parameters a, B/ 

X have quite different values in all runs. To judge their 
sampling distributions, histograms were drawn separately; 
Figures 6, 7 and 8 are the histograms of the estimated 
values of a, 3/ A, respectively. 

It is suggested by Figures 6 and 7 for A and 3/ 
that normality may be assumed with sample means 0.9246, 
0.05488 and sample variances 0.001065, 0.0002246, respec- 
tively. Their mean, median, trimean and midmean are pretty 
close to each other. But, in the case of a, the histogram 
(Fig. 8) is a quite different picture which looks something 
like an exponential distribution instead of normal distri- 
bution. The reason is clear, because the negative values 
of estimated a were not taken into account. They are 
physically infeasible. 

b . Nonlinear Least Squares Approach 

The least squares criterion used here can be 
stated formally and generally as follows : 

^ - 2 

minimize ^ (Y. - Y.) (4.10) 

i=l ^ ^ 

where N is the number of observations, Y^ is the fitted 
value of Y^ . In this case, expression (4.10) can be re- 
written as follows 

min y {Z.-- Z.(a,3,A)} (4.11) 

i=l ^ 
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TABLE I 



COMPUTER SIMULATION RESULTS FOR ESTIMATION PROCEDURE 
Sample size n = 50; Ntimber of Total Runs m= 400; 



Value of parameters a, B/ X Introduced: a=l, 3=0.05, X=0.1. 



Number of 
Runs 


X 


a 


6 


1 


0.0858 


2.2255 


0.0549 


2 


0.06591 


1.2745 


0.0779 


3 


0.08381 


9.498 


0.0688 


4 


0.0881 


0.9304 


0.0 39 


5 


0.1286 


3.7429 


0.0498 


6 


0.0883 


2.114 


0.069 


7 


0.0785 


2.8281 


0.0593 


8 


0.0945 


0.4073 


0.0584 


9 


0.089 


0.3833 


0.0452 


10 


0.1344 


1.3351 


0.0426 


11 


0.0677 


0.2111 


0.0729 


12 


0.1304 


2.4871 


0.0512 


13 


0.0714 


0.3211 


0.0623 


14 


0.069 


0.3547 


0.0657 


15 


0.0775 


0.4813 


0.0772 


16 


0.1618 


2.7416 


0.027 


17 


0.0827 


0.927 


0.0705 


18 


0.0816 


1.1663 


0.0836 


19 


0.0515 


0.3777 


0.0775 


20 


0.1089 


1.6439 


0.0522 


21 


0.0767 


0.8667 


0.0605 


22 


0 .0499 


0.3414 


0.0748 


23 


0.0569 


0.2118 


0.0608 


24 


0.0648 


0.2689 


0.0715 


25 


0:i415 


2.1178 


0.0518 


26 


0.1636 


1.3635 


0 .0186 


27 


0.0947 


6.101 


0.0664 
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TABLE I Cont. 



N\aitiber of 
Runs 


A 


a 


6 


28 


0.1063 


0.5235 


0.0623 


29 


0.0825 


2.6602 


0.059 


30 


0.13 


0.3412 


0.0303 


31 


0.817 


0.24 


0.0685 


32 


0.0843 


0.4377 


0 .0393 


33 


0 . 1445 


4.6941 


0.0442 


34 


0.0737 


0.2956 


0.0665 


35 


0.0 45 


0.1372 


0.0646 


36 


0.1108 


1. 1026 


0.0427 


37 


0.1277 


1.091 


0.052 


38 


0.0514 


0.1362 


0.062 


39 


0.1165 


0.6794 


0 .0469 


40 


0.0954 


3.7845 


0.0587 


41 


0.0563 


0.4123 


0.0738 


42 


0.0684 


0.2903 


0.0574 


43 


0.1263 


2.0 59 


0.0357 


44 


0.1335 


7.1595 


0.0226 


45 


0.0708 


0.465 


0.063 


46 


0.0872 


2.16 4 


0 .0657 


47 


0.0834 


0.3568 


0.0567 


48 


0.0429 


0.12 7 


0.0694 


49 


0.0428 


0.1295 


0.0639 


50 


0.0834 


0.6807 


0.0489 


51 


0.1141 


1.5349 


0.0566 


52 


0.1057 


1.1757 


0.0441 


53 


0.1104 


5.0692 


0.0433 


54 


0.1426 


0 . 8845 


0.0244 


55 


0.1084 


0.5775 


0 . 0505 


56 


0.1048 


1. 8399 


0.0562 


57 


0.0561 


0 . 3197 


0.0643 


58 


0.0607 


0 . 4564 


0.0607 


59 


0.0724 


0 . 3133 


0.0551 
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TABLE I Cont 



Nuinber of 
Rians 


A 


a 


6 


60 


0.0957 


0.9584 


0.0557 


61 


0.0804 


0.7378 


0.0423 


62 


0.0888 


0.1718 


0.0393 


63 


0.064 


0.3779 


0.9543 


64 


0.0324 


0.2202 


0.0791 


65 


0.0695 


0.2428 


0.0564 


66 


0.0903 


6.4194 


0.0542 


67 


0.0741 


0.4118 


0.0718 


68 


0.0462 


0.167 


0.0671 


69 


0 .0805 


1.0279 


0.0662 


70 


0 .0 469 


0.2157 


0.0693 


71 


0.0487 


0.2457 


0.0688 


72 


0.0682 


1.1847 


0.0665 


73 


0.0538 


0.1964 


0.0573 


74 


0.1516 


2.099 


0.0352 


75 


0.0763 


0.3067 


0.0679 


76 


0.0558 


1.5589 


0.079 


77 


0.0463 


0.2252 


0.0612 


78 


0.095 


5.6243 


0.0627 


79 


0.1903 


6.0 


0.0604 


80 


0.0822 


0.9989 


0.0738 


81 


0.066 


0.3422 


0.0926 


82 


0.0784 


0.3336 


0.0517 


83 


0 .0499 


0.1453 


0.062 


84 


0.0667 


0.1357 


0.064 


85 


0.0431 


0 .0 761 


0.0568 


86 


0.1485 


6.1674 


0.0262 


87 


0.1154 


1.0451 


0.0476 


88 


0.0043 


0.3493 


0 0648 


89 


0.0983 


1.3101 


0.0422 


90 


0.138 


6.7408 


0.0443 

— 
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TABLE I Cont 



Niimber of 
Runs 


A 


a 


6 


91 


0.1117 


0.7281 


0.058 


92 


0.065 


0.4284 


0.0482 


93 


0.0638 


1.2754 


0.0611 


94 


0.1603 


1.003 


0.0302 


95 


0.0436 


0.1647 


0.0 85 


96 


0.0905 


6.2455 


0.0703 


97 


0.1309 


1.4348 


0.0431 


98 


0.0787 


0.5245 


0.0686 


99 


0.1927 


7.801 


0.009 


100 


0.0683 


0.525 


0.0698 


101 


0.1163 


1. 7146 


0.0355 


102 


0.1005 


1.5795 


0.0463 


103 


0.1492 


1.5169 


0.0465 


10 4 


0.1511 


2.7243 


0.318 


105 


0.0823 


0 .1634 


0.0438 


106 


0.1493 


6.9467 


0.0437 


107 


0.1185 


9 . 4492 


0.0603 


10 8 


0.1039 


0 .589 


0.0406 


109 


0.0832 


0.3935 


0 .0612 


110 


0.059 


0.145 


0.0508 


111 


0.133 


8.6071 


0.0397 


112 


0.0877 


0.4751 


0.0448 


113 


0.0772 


0.2752 


0.0624 


114 


0.0659 


0.4782 


0.0583 


115 


0.0796 


0.2978 


0.0537 


116 


0.1128 


1.668 


0.0601 


117 


0.0974 


1.2931 


0.0515 


118 


0.1105 


0.3394 


0.0468 


119 


0.1216 


1.0529 


0.0655 


120 


0.1046 


2.0384 


0.0591 


121 


0.128 


3.2468 


0.0405 



45 



TABLE I Cont 



Number of 
Runs 


A 


a 


6 


122 


0.0818 


1.092 


0.0584 


12 3 


0.0924 


2.4295 


0.0439 


124 


0.0732 


0.1916 


0.0495 


125 


0.1284 


1.8515 


0.0407 


126 


0.1376 


2.6651 


0.0282 


12 7 


0.0877 


7.8429 


0.0703 


12 8 


0.0865 


1.1942 


0.0438 


129 


0 .1404 


2.8991 


0.0326 


130 


0.0452 


0.0978 


0.0663 


131 


0.0928 


3.0463 


0.0471 


132 


0.0668 


0 . 4602 


0.0551 


133 


0.2089 


0.7373 


0.0119 


134 


0.0946 


0.4011 


0.038 


135 


0.1352 


1.2162 


0.0284 


136 


0.0597 


0.2241 


0 .0539 


137 


0.106 


1.7821 


0.0354 


138 


0.0508 


0.2689 


0.0596 


139 


0.0522 


0.1244 


0.055 


140 


0.104 


0 . 8464 


0.0603 


141 


0.1002 


0.817 


0.0711 


142 


0.0692 


1.4214 


0.0744 


143 


0.0995 


0.5442 


0.0435 


144 


0.0697 


0.4966 


0.0487 


145 


0.0737 


0 .8995 


0.0708 


146 


0.059 


0.2998 


0.0654 


147 


0.1235 


1.1251 


0.045 


148 


0.1023 


0.615 


0.0349 


149 




0.9087 


0.0466 


150 


0.1242 


2.9267 


0.0356 


151 


0.0966 


0. 4989 


0.0616 


152 


0.0913 


2.5416 


0.0638 


153 


0.1301 


1.3897 


0.0434 
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TABLE I Cont 



Number of 
Runs 






6 


154 


0.1033 


2.0535 


0.0466 


155 


0.0638 


0.2219 


0.0594 


156 


0.1408 


3.6331 


0.0434 


15 7 


0.1311 


4.9526 


0.0495 


15 8 


0.0608 


0 . 469 


0.0462 


159 


0.0574 


0.5886 


0.0723 


160 


0.0998 


4.912 


0.0554 


161 


0.1114 


0.4898 


0.0363 


162 


0.0415 


0.1827 


0.0854 


163 


0.1586 


1.5257 


0.0397 


164 


0.0723 


0.9055 


0.0637 


165 


0.0984 


2.1521 


0.0675 


166 


0 . 10 11 


2.1535 


0.0573 


16 7 


0.093 


4.8488 


0.0555 


168 


0.1141 


4.2636 


0.0443 


169 


0.1953 


4.6512 


0.0452 


170 


0.0568 


0.8364 


0.07099 


171 


0.0927 


5.1303 


0.068 


172 


0.0632 


0 . 7399 


0.069 


173 


0.0536 


0.2732 


0.07138 


174 


0. 1033 


0.9515 


0.0787 


175 


0.0633 


1.6597 


0.0791 


176 


0.1414 


1.2547 


0.0318 


177 


0.1089 


4. 1786 


0.051 


178 


0.1165 


0.5441 


0.0432 


179 


0.1731 


4.2913 


0 . 0 19 2 


180 


0.0679 


0.1064 


0.0451 


181 


0.1682 


3.1417 


0.0332 


182 


0.0868 


0.6399 


0 .0685 


183 


0.0807 


1.046 


0.0607 


184 


0.0778 


1.185 


0.0694 
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I 



TABLE I Cont 



Number of 
Runs 


A 


a 


6 


185 


0.0769 


0.1999 


0.0009 


186 


0.0396 


0.1194 


0.064 


187 


0. 1012 


0.499 


0.485 


188 


0.1389 


9.70 85 


0.0334 


189 


0.0665 


1.1661 


0.0599 


190 


0.0471 


0.1374 


0.0643 


191 


0.1316 


4.1676 


0.0259 


192 


0.0687 


0.6692 


0.0928 


19 3 


0.0866 


0.5419 


0.0484 


194 


0.0966 


1.5 79 


0.0524 


195 


0.0773 


1.0102 


0.0668 


196 


0 . 1076 


0.7413 


0.0704 


197 


0.0857 


0.8715 


0.0548 


198 


0.0803 


1.110 


0.0549 


199 


0.052 


0.1396 


0.0598 


200 


0.0757 


0.9558 


0.0622 
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FIG. 6. Histogram for 
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FIG. 7. Histogram for 
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Fig. 8. Histogram for 



The goal is to choose values a, 3 and X which minimize 
the expression (4.11) . In the linear case, this can be done 
using elementary calculus by taking the partial derivatives 
with respect to a, 3 and X , setting each of them equal 
to zero and solving the resulting linear equations (Normal 
equations) . But the present problem is different because 
the expression (4.11) is highly nonlinear in the parameters, 
and indefinite in terms of convexity. This is why a non- 
linear programming technique is used. 

There are a few general approaches to the solu- 
tion of the nonlinear estimation problem. One of them is 
the direct optimization approach. 

Specifically, equation (4.9) is considered as 



follows : 



where 



P1P2 



(p) 



~ 11 + p2£(p)>ii + p3£(p)y 



P 



1 



1 

X 



P2 



a 

X 



P3 



3 

X 



(4.12) 



In equation (4.12), z(p) and e (p) are known values, 
are unknown variables. If it is rewritten as the sum of 
squared errors the objective function will be ; 
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Now nonlinear optimization problems can be stated such that 



N ( 

Min I .J z (P. ) 
i=l ( ^ 



PlP2£^(Pi) 



(1 +P2£(P^)) (1 + P3S(P^)) 



Subject to 



Pl^O, P2 - ^ ' P32.0 



In order to perform the optimization, i.e. to 
solve this nonlinear optimization problem, the GRG package 
(Generalized Reduced Gradient) [6] was used; it is very 
convenient for this case. Ten runs were made with five 
different initial points, both with analytically computed 
derivatives and with numerically computed derivatives. The 
first four initial points were selected arbitrarily. The 
last one was picked such that: 



Pi 




P2 



P 



3 



^(1) 

X 

^(50) 

✓N 

X 



(4.14) 



where ^(j_) smallest number in the sample, 

is the largest number in the sample, and X is the natural 
logarithm of p = 0.5. divided by the median of the sample. 
Initial points are shown in Table II. 

The optimization results of a, and p 

associated with Pj^, P 2 and p^ are tabulated in 
Table III and Table IV, respectively. 
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As it is seen from Tables III and IV, the program 
was terminated at different optimal points. Because the 
objective function is highly nonlinear and apparently in- 
definite (neither convex nor concave) , it may sometimes 
have stopped at a local minimum rather than at a global 
optimum point. But if the results are compared to the 
3 -percentile results, the nonlinear least square estimates 
show much greater accuracy and seem more consistent. 



TABLE II 

INITIAL GUESS POINTS 



# of 
Points 


Pi 


P2 


P3 


I 


50.0 


2.0 


30.0 


II 


10.0 


10.0 


0.5 


III 


80 .0 


80.0 


5.0 


IV 


1.0 


1.0 


1.0 


V 


Accordin( 


g to the e< 


5 . (4.14) 
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Computed by using 



TABLE III Cont. 
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TABLE IV 



Actual values A = 0.1, a = 1.0, 3 = 0.05 



# of 
Runs 


Initial* 


A 


a 


B 


Point 


A 


B 


A 


B 


A 


B 




I 


0.1202 


0.1203 


343.512 


676. 800 


0.0400 


0.0400 




II 


0.1199 


0.1203 


129.591 


968.150 


0.0401 


0.0401 


1 


III 


0.1168 


0.1203 


8.898 


1506.926 


0.0409 


0.0400 




IV 


0.1158 


0.1204 


6.577 


874.957 


0.0412 


0.0385 




V 


0.1203 


0.1201 


4.182 


955.629 


0.0400 


0.040,1 




I 


0.0984 


0.0863 


7.534 


1.297 


0.0478 


0.0510 




II 


0.0865 


0.0863 


1.306 


1.296 


0.0509 


0.0510 


2 


III 


0.0977 


0.0863 


6.158 


1.297 


0.0480 


0.0510 




IV 


0.0798 


0.0863 


0.810 


1.292 


0.0527 


0.0511 




V 


0.0500 


0.0863 


0.154 


1.296 


0.0533 


0.0511 




I 


0.1390 


0.1434 


8.123 


392.550 


0.0343 


0.0331 




II 


0.1430 


0 . 1434 


138.531 


1520.067 


0.0332 


0.0332 


3 


III 


0.1440 


0.1436 


204.078 


1465.586 


0.0328 


0.0330 




IV 


0.1434 


0.1433 


414.903 


1177.650 


0.0331 


0.0 332 




V 


0.0438 


0.1434 


0.055 


2056 .280 


0.0515 


0.0331 




I 


0.0878 


0.0843 


0.118 


0.10 7 


0.0394 


0.0397 




II 


0 .0890 


0.0873 


0. 12 3 


0.116 


0.0392 


0.0393 


4 


III 


0 .0891 


0.0856 


0.12 3 


0.111 


0.0393 


0.0395 




IV 


0.0870 


0.0829 


0.115 


0.10 3 


0.0394 


0.0398 




V 


0.0899 


0.0860 


0.125 


0.112 


0.0390 


0 . 0 39 4 




I 


0.0507 


0.0508 


4.944 


0.509 


0.0530 


0.0775 




II 


0.0695 


0.0509 


0.847 


0.513 


0.0988 


0.0775 


5 


III 


0.0692 


0.0508 


3.882 


0.508 


0.0749 


0.0775 




IV 


0.0540 


0.0508 


0.826 


0.509 


0.0717 


0.9775 




V 


0.0676 


0.0507 


0.6 72 


0.506 


0.1034 


0.0776 
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TABLE IV Cont 



# of 
Runs 


Initial 

Point 


A 


a 




A 


B 


A 


B 


A 


B 




I 


0.1384 


0.0994 


9.0 72 


0.570 


0.0494 


0.0591 




II 


0.1019 


0.0996 


0.628 


0 .57 6 


0.0585 


0.0590 


6 


III 


0.1363 


0.0998 


6.105 


0.5 79 


0.0499 


0.0589 




IV 


0.1079 


0.0998 


0.810 


0.5 82 


0.0571 


0.0590 




V 


0.1272 


0.0994 


2.439 


0.571 


0.0525 


0.0590 




I 


0.1158 


0.0738 


6.780 


0. 301 


0.0405 


0.0502 




II 


0.0755 


0.0738 


0.324 


0.301 


0.0496 


0.0502 


7 


III 


0.0752 


0.0743 


0.319 


0.309 


0.0499 


0.0501 




IV 


0.0757 


0 .0739 


0.328 


0. 30 8 


0.0499 


0.0501 




V 


0.1188 


0.0738 


28.839 


0.301 


0.0397 


0.0502 




I 


0.0753 


0.0748 


7.325 


5.783 


0.0645 


0.0631 




II 


0.0739 


0.0740 


4.647 


4.76 3 


0.0634 


0.0634 


8 


III 


0.0735 


0.0735 


4.219 


4.221 


0.0635 


0.0635 




IV 


0.0737 


0.0741 


4.361 


4. 831 


0.0634 


0.0633 




V 


0.0170 


0.0010 


0.033 


0.0008 


0.0724 


0.0634 




I 


0.1257 


0 . 1252 


11.646 


9. 302 


0.0275 


0.0276 




II 


0.1252 


0. 1253 


10.181 


9.374 


0.0275 


0.0275 


9 


III 


0.1240 


0.1246 


6.983 


8.598 


0.0278 


0.0278 




IV 


0.1254 


0.1254 


11.010 


10.296 


0.0276 


0.0275 




V 


0.1017 


0.1255 


0.754 


10.407 


0.0330 


0.0275 




I 


0.1387 


0.1438 


5.688 


25 . 170 


0.0326 


0.0313 




II 


0.1441 


0.1428 


28.830 


29.290 


0.0312 


0.0310 


10 


III 


0.1400 


0.1440 


6.290 


27.747 


0.0326 


0.0313 




IV 


0.1437 


0.1441 


21. 370 


26.590 


0.0313 


0.0312 




V 


0.1415 


0.1427 


10.030 


14.610 


0.0310 


0.0316 
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C. TEST PROCEDURE FOR PARAMETERS 

In order to determine whether the model reasonably fits 
the data a simple test procedure is applied. It is to 
exhibit the simple comparison of the actual points z(P) 
and the estimated values z(P), using the estimated 
parameters . 

Basically z (P) can be obtained as follows; 

x(p) = - i S,n(l-p) (4.15) 

X 

Then if it is substituted in expression (2.37), z (P) will 
be 



z(P) 



/V /\ O 

ax (P) 

(1 + ax(P) ) (1 + Bx(P) ) 



(4.16) 



In equation (4.16) , all variables are known so estimated 
quantiles can be easily obtained. 

In Fig. 9 actual values and estimated values are plotted. 
Also Fig. 10 displays residuals, which are the difference 
between actual and estimated values. 
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Fig, 9. Comparison of Actual Values and Estimated Values 



o 

(O 
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Fig. 10. Residuals of Fitting 



V. NUMERICAL APPLICATIONS TO TWO SETS 
OF REAL DATA 

In this section the two failure data sets studied come 
from these sources: Oral irrigator [7] and human life [8]. 

They are used for numerical examples in that these data are 
fitted using the model (2.6). 

A. ORAL IRRIGATOR 

The data used was obtained from the Commun. Statist. - 
Theor. meth. , Colvert and Bordman [7] . The data was col- 
lected such that 100 oral irrigators were placed on the test 
was terminated at 700 time units. During the life test, 98 
oral irrigators failed; another 2 oral irrigators survived. 
The ordered observed times to failure of the oral irrigators 
is tabulated in Table V. 

Using the data of Table V, the parameters a, 8, and X 
are estimated by means of the 3-percentile approach. To 
demonstrate the differences between estimated values, 
various different combinations of p values were used. The 
results are shown in Table VI. 

Examination of Table VI indicates that the 3-percentile 
approach may produce estimates having great differences for 
different percentile values, except here in the case of the 
first two combinations. Also the first two estimations seem 
to have accepatable limiting age. Moreoever, the last two 
combinations do not fit well, as judged from the residuals. 
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TABLE V 



TIME TO FAILURE OF ORAL IRRIGATORS 



1. 75 


7.02 


7.58 


9.76 


15.02 


15.5 7 


17. 39 


19.55 


22.47 


23.24 


23.96 


25.05 


32.44 


36.87 


42.76 


43.14 


46.95 


56.33 


58.99 


59.08 


60.37 


61.01 


• 77.86 


86.45 


88.50 


103.06 


104.34 


105. 85 


117.46 


120.11 


122.28 


122.61 


129.31 


130.42 


137.57 


142.27 


142.98 


148.29 


150. 79 


151.21 


155.62 


157.93 


160. 72 


169. 79 


186.26 


197.60 


224. 83 


233.64 


242.07 


256.86 


260.77 


261.68 


277.99 


283.95 


288.94 


295.48 


314. 76 


316.06 


332.07 


339.46 


362.61 


369.47 


370. 74 


491.06 


403. 39 


414.78 


426.71 


459.62 


455.84 


457.94 


466.61 


468.64 


469.09 


476.42 


481.41 


481.82 


488.15 


490.06 


493.67 


494.38 


503. 72 


508.93 


509.01 


418.32 


532.29 


534.62 


545.23 


547.41 


558.41 


571.10 


585.52 


589.11 


592.93 


607.15 


623.15 


647.91 



TABLE VI 



ESTIMATED VALUES OF a, 8, A BY USING 3-PERCENTILE 



Pi 


P 2 


P 3 


A 


/s 

a 


3 


1/3 


.1 


.5 


.9 


0.00153 


0.00634 


0.001028 


972.76 


.25 


.5 


.9 


0.00158 


0.00710 


0.001018 


982.31 


.1 


.5 


. 75 


0.00239 


0.01835 


0.000228 


4385.96 


.25 


.5 


.75 


0.00273 


0.04514 


0.000077 


12987.01 
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Next the nonlinear optimization approach was applied to 
estimate parameters a, 3 and X associated with 
snd P 3 ' hoth using the derivative and without derivative by 
taking into account 98 ordered values. The results >are 
tabulated in Tables VII and VIII for p^, p^, P 3 and X 
a, 3 respectively. Some initial points in Table II are 
used. Table VII indicates that the estimates obtained from 
the nonlinear estimation method are more consistent than 
those from the 3-percentile approach. Also, the boundary 
points, 1/3, seem reasonable. Finally, comparison of the 
actual values of z(p) and estimated values of z(p) indi- 
cate that the fitting is reasonable. The residuals of 
fitted model and the comparison of the actual and estimated 
values are plotted in Fig. 11 and Fig. 12. 

B. HUMAN LIFE (MORTALITY) DATA 

The data used in this example was obtained from the 
1969-71 life table [ 8 ] for white females in the United States 
The table has been prepared from a history of 100,000 persons 
the number of surviving, and the number dying has been given 
for each age interval. Look at Table IX. 

1. Estimation by Using 3-Percentile Approach 

The life data is first fitted to the model (2.37) 
by using the 3-percentile approach. Certain pth percentile 
values are selected and used in the estimation process. 

These results are tabulated in Table X. Investigation of 
the last eight combinations of percentiles in these tables 
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ESTIMATED PARAMETERS FOR ORAL IRRIGATORS 
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Fig, 11. Comparison of Actual Values and Estimated Values 
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Fig. 12. Residuals of Fitting 



indicates that the limiting age, 1/B , is approximately 110, 
which is reasonable. However, the first combination of 
percentiles gives a good fit in terms of the sum of the 
squared errors. But for these limiting age is reduced to 
about 95 years. This apparently means that the fit of the 
model does not represent the age above 95. Thus it can 
be said that there is a trade-off between 1/3 and the sum 
of squared errors. The present model simply does not seem 
to fit the mortality data very well. 

2 . Estimation By Using Nonlinear Least Squares, and 
Application of a Nonlinear Programming Algorithm 
In the previous section we discussed the fact that 
one of the problems encountered in the 3-percentile estima- 
tion procedure was the high uncertainty of fitting, as 
measure by the sum of squared errors. In order to reduce 
this variability, a nonlinear estimation process is applied. 
Previously it was noted that the GRG package can be used 
with analytically computed derivatives and without 
derivatives. If it is used with analytically computed 
derivatives, it is necessary to provide another subprogram 
by the user. Otherwise the GRG package will provide the 
derivative, computed numerically and automatically. 

The derivative of the objective function (4.13) is 
the normal equations such that 
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TABLE IX 

MORTALITY DATA FOR WHITE FEMALES, 1969-71 



Time 

Interval 


# of 
Survivors 


# of 
Dying 


Time 

Interval 


# of 

Survivors 


# of 
Dying 


0- 1 


100,000 


15 32 


27-28 


97,165 


70 


1- 2 


99,468 


100 


28-29 


97,095 


73 


2- 3 


98,368 


65 


29-30 


97,022 


77 


3- 4 


98,303 


54 


30-31 


96,945 


81 


4- 5 


98,249 


46 


31-32 


96 , 864 


87 


5- 6 


98,203 


39 


32-33 


96,777 


93 


6- 7 


98,164 


35 


33-34 


96,684 


101 


7- 8 


98, 129 


32 


34-35 


96,583 


109 


8- 9 


98,097 


29 


35-36 


96,474 


118 


9-10 


98,068 


26 


36-37 


96,356 


12 8 


10-11 


98,042 


24 


37-38 


96,228 


141 


11-12 


98,018 


24 


38-39 


96,087 


153 


12-13 


97,994 


25 


39-40 


95,932 


170 


13-14 


97,969 


30 


40-41 


95,762 


185 


14-15 


97,939 


37 


41-42 


95,577 


201 


15-16 


97,902 


45 


42-43 


95,376 


220 


16-17 


97,857 


54 


43-44 


95,156 


242 


17-18 


97,803 


60 


44-45 


94,914 


263 


18-19 


97, 743 


62 


45-46 


94,649 


291 


19-20 


97,681 


63 


46-47 


94,358 


317 


20-21 


97, 618 


63 


47-48 


94,041 


344 


21-22 


97,555 


63 


48-49 


93,697 


372 


22-23 


97,492 


63 


49-50 


93, 325 


401 


23-24 


97,429 


64 


50-51 


92,924 


433 


24-25 


97,365 


66 


51-52 


92,491 


469 


25-26 


97,299 


66 


52-53 


92,022 


506 


26-27 


97,233 


68 


53-54 


91,516 


546 
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TABLE IX Cont 



Time 

Interval 


# of 

Survivors 


# of 
Dying 


Time 

Interval 


# of 

Survivors 


# of 
Dying 


54-55 


90,970 


587 


82- 83 


41,215 


3,586 


55-56 


90,383 


6 32 


83- 84 


32,629 


3,589 


56-57 


89,751 


680 


84- 85 


34,040 


3,550 


57-58 


89,071 


730 


85- 86 


30,490 


3,495 . 


58-59 


88,341 


781 


86- 87 


26,995 


3,425 


59-60 


87,560 


834 


87- 88 


23,570 


3,2 86 


60-61 


86 , 726 


911 


88- 89 


20,284 


3,0 72 


61-62 


85,835 


953 


89- 90 


17,212 


2,806 


62-63 


84,882 


1,022 


90- 91 


14,406 


2,531 


6 3-6 4 


83,860 


1,098 


91- 92 


11,875 


2,262 


6 4-6 5 


82,762 


1,183 


92- 93 


9,613 


1,982 


65-66 


81,579 


1,275 


93- 94 


7,631 


1,694 


66-67 


80 ,30 4 


1,375 


94- 95 


5,937 


1,411 


67-68 


78,929 


1, 486 


95- 96 


4,526 


1,145 


68-69 


77,443 


1,607 


96- 97 


3,381 


905 


69-70 


75,836 


1,735 


97- 98 


2,476 


696 


70-71 


74,101 


1,862 


98- 99 


1,780 


524 


j 71-72 


72,239 


1,993 


99-100 


1,256 


384 


72-73 


70 ,246 


2,141 


100-101 


872 


277 


73-74 


68,105 


2,312 


101-102 


595 


195 


74-75 


65,793 


2,503 


102-103 


400 


135 


75-76 


63,290 


2,693 


103-104 


265 


92 


76-77 


60,597 


2,872 


104-105 


173 


61 


77-78 


57,725 


3,0 39 


105-106 


112 


41 


78-79 


54,686 


3,187 


106-107 


71 


26 


79-80 


51,499 


3,317 


10 7-10 8 


45 


17 


80-81 


48,182 


3,434 


108-109 


28 


11 


81-82 


44,748 


3,533 

— t- 


109-110 


17 


6 
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TABLE X 



ESTIMATED VALUES OF a, 3, A FOR MORTALITY DATA USING 

3 -PERCENTILE APPROACH 



Percentile 

Values 


/\ 

A 


ys 

a 


/N 

3 


/s 

-.1/3 


22 


56 


92 


0.000674 


0.1190 


0.010550 


94.7 


20 


56 


92 


0.000550 


0.0404 


0 .010520 


95.0 


17 


56 


92 


0.000508 


0.0302 


0.010550 


94.7 


16 


56 


92 


0.000440 


0.0199 


0.010550 


94. 7 


17 


56 


109 


0.000862 


0.1676 


0.009069 


110.2 


18 


56 


109 


0.000886 


0.2389 


0.009068 


110.2 


19 


56 


109 


0.000908 


0.3708 


0.009067 


110.2 


17 


56 


110 


0.000876 


0.1805 


0.008991 


111.2 


18 


56 


110 


0.000899 


9.2613 


0.008989 


111.2 


19 


56 


110 


0.000920 


0.4188 


0.008988 


111.2 


20 


56 


110 


0.000941 


0.9070 


0.008988 


111.2 
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IS 

9p 



1 ' ^ p2£(Pi)) (1 + P3£(Pi)) I 



N 

2 I 



P2^ (P^) 



Plp2£ (P-) ) 

X < 2 ( P ) — ^ i ' 

i (1 + P2£(P£))(1 + (Pj^) ) I 



(4.17) 



as 

3P2 



_ 2 I jpi£^(Pj^)(l + p2£(P^)) - £ (P^) pj^p2£^(P^) 

i=l j (1 + P2£(P^))^ (1 + P3£(P^)) 



I ^1^2^ 

( ^ (1 + P2£ (P“)) Tl + P,£ (P,)) 



(4.18) 



3 '"i' 



as 

ap 3 



+ 2 



I 

i=l I (1 + P2£(P.))"^ (1 + P3E(P.)) 



z(p^) 



Pi p 



1^2 



£^ (P, ) 



(1 + P^£ (P,) ) (1 + 



P2£ ) U + P3£ (P^) ) j 



(4.19) 



N is 110 in the case of human mortality data. 

The results of estimation were shown in Tables XI 
and XII for p^ and a, 8f A respectively. Also Fig. 13 
and Fig. 14 demonstrate the actual and estimated values and 
residuals of fitting. 

Examination of Table XI and Table XII indicates 
that the objective function values (the sum of the squared 
errors) of the fitting are smaller than the objection func- 
tion values of 3-percentile approach. Also the estimated 
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parameter values 3 and X associated with and 

are very close to each other. But estimated values of a 
show considerable difference for the various initial values. 
However, there are two important facts to notice. First, 
when ot increases significantly, the objective function 
values remain almost constant. That is, it is not effected 
significantly on the objective function values in this case. 
Second, this analytical model (2.6) does not represent 
deaths beyond the age 95. However, modifications can be 
made in the model (2.6) for this kind of difficulty. Some 
ideas will be discussed in a later section. 



3 . Model Modifications 

In the previous section it was noted that the model 
(2.6) has relatively great errors and does not accurately 
describe probability of death at age greater than 95 in the 
hiaman life example. This means that the hazard fvmction 
h (z) increases too rapidly in the wearout period. If it 
is possible to slow the rate of increase of hazard perhaps a 
better result can be obtained. 

In Section II, it was stated that the function 
R(x) describes the wearout period in bath tiob type curves; 



R(x) 



1 

1 + 3x ' 



3 > 0 



If a new parameter, y, is defined which is between 

th 

0 and 1 and R' (x) is now defined to be the y power of 
R(x) : 
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R' (X) 



1 



8>0, 0<Y<1 



(4.20) 



(1 + 8x)^' 

then R' (x) will provide for a slowly increasing, rather 
than rapidly increasing, in wearout period. With this 
revision the model (2.6) now becomes 



or 



z 



2 



G(x) = xL(x) • R' (x) 
2 

OJC 

(1 + ax) (1 + 8x) ^ 



(4.21) 



(4.22) 



where a > 0, 6 > 0 and 0 < y < 1- 

However, after we put another variable in the model, 
it will be too difficult to handle in the previous estimation 
technique for obtaining values of the parameters a, 6/ X 
and Y because of nonlinearity and indefiniteness of the 
expression (4.22). 

For this reason, y will be assumed constant in the 
previous estimation procedure. Then it will be computa- 
tionally convenient. Actually, the estimation procedure 
does not require any change. The nonlinear least square 
estimation approach will simply be used for different values 
of y. Table XIII and Fig. 15 demonstrate the parameters 
estimated and the resulting objective function values. From 
Table XIII, it can be easily seen that the new result has a 
smaller sum of squared errors. The best value of y seems 
to be near 0.95. 
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ESTIMATED PARAMETERS FOR HUMAN LIFE BY USING SAME INITIAL POINTS 
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Fig. 13. Comparison of Actual Values and Estimated Values 
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Fig. 14. Residuals of Pitting 
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Fig. 15. Plot of Different y Values vs Objective Values 



TABLE XIII 



ESTIMATION VALUES OF a, 6, A 
USING CONSTANT y 



# of 
Run 


Y 


X 


a 


6 


Objective 

Value 


1 


1.000 


0.000834 


0.42950 


0.01029 


4,335.89 


2 


0.97 


0.000512 


0.04051 


0.01182 


3,119.85 


3 


0.95 


0.000363 


0.01928 


0.01314 


2,494. 10 


4 


0.93 


0.000740 


1.06371 


0.01352 


2,999.11 


5 


0.90 


0.000685 


1.81283 


0.01554 


3,267.47 


6 


0.85 


0.000592 


5.74169 


0.02031 


5,099.37 


7 


0.7 


0.000319 


3.17467 


0.04580 


15,310.06 
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C. CONCLUSIONS 



Simple analytical hazard models have been developed 
and fitted to situations (data) that exhibit bath tub 
shaped hazard functions. That is, failure rates may be 
high at early ages ("infant mortality") , constant at 
intermediate ages, and high again for later ages ("wearout") . 
The procedure emphasizes representations of the inverse 
distribution function; simulation is thus facilitated. 

The failure time distributions so derived should be 
useful in analyzing maintenance and replacement policies. 

A least squares technique for fitting the hazard models 
to data are suggested and applied. 
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